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FOREWORD 


This  paper  was  prepared  by  Mr.  J.  B.  Cheek,  Jr.,  Dr.  N.  Radhakrishnan, 
and  Mr.  F.  T.  Tracy  of  the  Computer  Analysis  Branch,  Automatic  Data  Proc- 
essing Division,  U.  3.  Army  Engineer  Waterways  Experiment  Station  (WES), 
for  presentation  to  the  1971  Army  Numerical  Conference  sponsored  by  the 
Army  Mathematics  Steering  Committee  and.  hosted  by  the  Department  of  Defense 
Computer  Institute.  The  paper  was  reviewed  and  approved  for  presentation 
and  publication  by  the  Office,  Chief  of  Engineers,  U.  S.  Army. 

The  work  was  conducted  during  the  period  Dec  19&7  to  Apr  1971  under 
the  general  supervision  of  Mr.  D.  L.  Neumann,  Division  Chief.  It  is  based 
n work  done  for  several  of  the  WES  technical  divisions  and  the  Lower  Missis- 
sippi Valley  Division  in  connection  with  a number  of  engineering  projects. 

During  the  period  in  which  this  paper  was  prepared,  COL  Ernest  D. 
Peixotto,  CF , was  Director  of  WES;  Mr.  F.  R.  Brown  was  Technical  Director. 
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SUMMARY 


This  paper  was  prepared  to  familiarize  practicing  scientists  and  en- 
■ ■ rs  with  the  cubic  spline  interpolation  technique  as  a possible  tool  in 
:rv<  fitting  for  computer  programs  for  which  more  commonly  used  techniques 
-ay  be  unsuitable  or  of  limited  value.  The  spline  technique  is  compared 
v more  • mm  on  methods,  specifically  piece-wise  linear  and  polynomial, 

■md  '.xamjles  of  applications  of  the  technique  to  engineering  problems  are 
presented.  Appendix  A contains  the  mathematical  derivation  of  the  equa- 
*'.ons  iefining  the  spline  function,  and  Appendix  B contains  a compact 
RTRA\  fitting  and  interpolating  program. 

The  interpolating  spline  curve- fitting  technique  has  three  primary 
advantages : 

a.  The  spline  passes  through  all  data  points. 

b.  The  first  and  second  derivatives  of  the  spline  are  continuous 
at  all  px)ints . 

c.  The  spline  can  be  easily  modified  to  satisfy  new  or  addi- 
tional data. 

he  experience  of  the  Waterways  Experiment  Station  (WES)  in  applying  spline 
fec.Miiques  to  engineering  problems  has  indicated  that  these  advantages  out- 
weigh the  additional  storage  and/or  computation  time  requirements  of  the 
technique  in  -any  applications. 

Since  the  spline  function  is  required  to  pass  through  all  data 
points,  erratic  derivative  behavior  may  result  from  experimental  error 
when  the  data  points  are  numerous  and  closely  spaced.  Trial  and  error 
.methods  for  smoothing  such  functions  exist,  but  they  are  time  consuming. 

WES  experience  has  indicated  that  acceptable  results  can  generally  be 
obtained  by  simply  selecting  a more  limited,  more  widely  spaced  set  of  the 
lata  points  to  which  to  fit  the  curve. 
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pony  computer  pr  grams  apj  .icd  • engineering  research  and  design 
list  i'  . th<  s rs  teristics  f nonlinear  physical  systems.  The 
line  t'fers  utstaniing  advantages  >ver  the  interp  lating  methods 
. used  in  this  class  of  problems.  It  is,  therefore,  the  purpose  of 
• is  taper  t point  ut  the  shortcomings  f several  methods  and  show  how 
split.1  techniques  may  be  used  t advantage.  This  purpose  is  approached 
■ ,r  ugh  iis  -ussion  and  examples  in  language  and  subject  that  are  meaningful 
research  and  design  engineer  in  order  to  bring  to  the  practicing 
scientist,  and  engineer  an  assurance  that  the  cubic  spline  formulation 
: •'ers  a powerful,  practical  modeling  find  interpolating  method  for  use  in 
is  • mputer  codes. 


Commonly  Used  t'urvo-F i t.t  i ng  Pechnj  iuei 


2.  Numerical  techniques  and  digital  computers  are  being  applied  to 

nuiber  f :ivi3  engineering  problems.  This  is  principally 
:.  1 th<  fli  - i:  i i i ■ iff  rd<  i by  the  numeri cal  procedures,  in  that  the 

design  and  research  engineer  can  easily  specify  complicated  boundary  con- 
ditions and  use  nonlinear  material  properties.  The  finite  dil t erence  and 
finite  element  methods  are  excellent  examples  of  this  growing  application 
area.  The  valid  description  (modeling)  >f  the  nonlinear  properties  to  the 
omputer  program  is  a primary  consideration  and  is  the  major  concern  ol 
this  paper. 

3.  Problems  of  this  type  require  that  nonlinear  relationships  be- 
tween physical  quantities  be  available  to  the  solution  process  in  either 
functional  or  analytical  form.  This  is  necessary  to  answer  questions  like 


r.  a s • ra in . wha t is  th<  tress ; r given  water  eleval 
what  is  the  flow  (discharge).  Such  rolat  ionshij.s  are  normally  avai  v 
only  as  lata  joints,  it  i.:  therefore  necessary  to  use  ram'-  kind  of  -ur.  - 
:'i  1 1 i:.p  technique  in  the  *omputer  program  t represent  the  phy.ii  -a 
tv-t  only  at  'he  data  points,  hut  in  the  interval.'  t-  tween  lata  T-c 

U.  Th*  r-  are  many  -'urve-f  itting  methods,  t u*  • 3 • no  id  a : 

'■ons  quently.  a major  difficui  ' j in  level oping  the  solution  pr  e« 
s looting  one  •■urve-f  i ♦ ? ing  technique  that,  is  l *>st.  y.u  i ted  to  th-  :r 
hand . 

t.  Two  of  the  most,  commonly  used  curve-fitting  techniques  are  . 
v;ise  linear  and  polynomial  methods,  while  hyperiolic  fune'ion  and  ■ 
special  purpose  representations  nr--  occas  tonal  1 y us«-d.  Th-  r<»  ar  . • : . 

s rious  disadvantages  in  using  the  linear  and  polynomial  methods, 
limitations  are  presented  in  the  following  paragraphs  to  help  the  rea-.i  r 
appreciate  the  advantages  of  'he  spline  method. 

Piece-wise  linear  interpolation 

’ • disadvantage;; . Given  a series  of  data  points  that  , for  --ncd. 
ing  purposes,  ’exactly"  represent  a physical  situation,  there  is  .cron/ 
motivation  to  use  piece-vise  linear  interpolation  between  p-Lnt  pairs. 

. his  is  especially  tempting  when  additional  data  points  are  easy  * o a ■ : 

because  the  nonlinearity  can  be  modeled  (to  the  extent  requir'd ) sin-/. 

pecifying  additional  points  f r the  nonlinear  ("curvy"  regions.  It 
appears  that  the  only  disadvantage  is  th°  additional  computer  men:  -r\  r-  - 
quired  for  the  points. 

7-  The  serious  objection  of  discontinuous  derivatives  develo:  wi.- 

modeling  observed  physical  behavior  and  calculating  rates  of  change  fder'.v-- 
tives)  from  the  piece-wise  linear  model.  This  difficulty  caused  by  dis  ■ 
tinuous  derivatives  is  illustrated  in  the  following  example. 

8.  The  nonlinear  properties  of  soil  in  a finite  element  method  (FEM 
olution  may  be  represented  with  a set  of  stress-strain  points,  as  shown 
’■'  *•  solution  process  requires  interpolation  of  a stress  value 

for  any  strain  value.  Also  required  is  the  modulus  of  elasticity  for  those 
same  strain  values  (the  modulus  being  dY/dX  , the  slope  of  the  stress- 
rain  curve).  Fig.  l also  includes  the  plot  of  modulus  versus  strain. 
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h<  abrupt  -liangt  in  modulus  a 
from  regioi  1 I regioj  . 


KF.M  s o 1 u t ion  proc" 


•,r<'e  ■•  arao ter i;: t. i -s  of  a 'epical  non- 
M pi  : th<  ution  pr  lui 


solving  a series  of  incremental]  1 • >ad- 
ens,  • s 1 during  the  load  cycle,  the 
iu'us  of  • ach  soil  element  is  dependent  on 
strain  a*  ’hat  element,  and  (c)  the  strains 


A/iA  STRAIN 


Fig.  1.  Piece-wise 
linear  stress-strain 
approximation 


1:.  :r-as'  incrementally  from  an  ini’ial  value  as 
aiiitional  load  cycles  are  taken. 

10.  Thus , t hose  elements  having  any  value  of  strain  in  region  A will 
save  ore  modulus  valu<-,  while  those  having  any  strain  value  in  region  B 
■xhibit  a different  modulus.  Consequently,  as  the  element  strain  progresses 
from  region  A to  F,  the  solution  process  sees  an  abrupt  change  in  ‘hat  ele- 
ment’s modulus.  Such  behavior  is  n t.  -haracterietic  of  soil  modulus  versus 
train;  L.e. , the  mod--.;  for  modulus  versus  strain  is  obviously  invalid. 

The  effect  of  this  abru]  t changt  1 instal  ility  in  the  numerical  process 
and  questionable  resul’s. 

Linear  interpolation  is  an  -’asy  method  t.o  implement  in 
the  program.  It  can  accurately  model  the  observed  behavior,  provided  a 
sufficient  number  of  data  points  are  available;  but  it  fails  completely  in 
model  in.’  derivatives  and  is  wasteful  of  computer  memory  when  stringent 
limits  art-  placed  on  allowable  errors. 
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dynomia]  interpolation 


12.  Disadvantages . Because  linear  interpolation  has  discontinuous 
first  derivatives,  researchers  have  turned  to  higher  (N  ) degree  poly- 
nomials which  have  the  much  desired  continuity  of  their  first  through 
(N  - l)  ^ derivatives.  In  so  doing  they  discovered  difficulties  and  de- 
ficiencies associated  with  polynomial  interpolation,  some  of  which  are  out 
1 ined  bel ow : 

a.  Polynomials  of  degree  N do  not  always  pass  through  every 
data  point  (if  there  are  more  than  N + 1 data  points). 
This  is  objectionable  when  the  observed  phenomenon  is  such 


• method  ( FEM 
s , as  shown 
stress  value 
■ity  for  those 
he  stress- 


that,  -la*  a points  can  he  considered  "exact"  (i.<  la*  a 
points  whose  protable  error  of  position  is  small  ar.  i *• 
polynomia’  fit  does  not  mee*  *his  error  crition). 

b.  Oscillations,  as.  many  as  (N  - l)/2  , may  occur  betw 
the  firs*  and  las*  daia  j»ijit . ’ueh  oscillati  : ■■  ■■. 

outside  the  rang*  of  the  known  behavior  of  *h  . * •• 
modeled  or  may  incorrectly  model  derivative  behavior. 

It  is  extremely  difficult  t predict  the  • -ye  rail  ■ • a ■■ 
an  Nth  degree  polynomial  when  a data  point  is  added, 
leted,  or  "adjusted"  (again  assuming  then  ar  me-  • • a 
dat.a  points).  This  difficulty  ’ransformr  ;urve-:'l*- 

process  from  a science  to  an  ar* . 

d.  Although  derivatives  of  polynomials  are  eordinuous. 

not  assume  that  they  are  representative  of  the  i iea. 
situation  being  modeled. 


13-  Summary.  Polynomial  interpolation  methods  can  l"  su-’c-ss-'. 
applied  to  many  nonlinear  problems,  but  a high  level  of  edu^at  i n,  •<’. 
and  experience  is  required  of  the  curve-fitting  prac'itior:  r (the  - qua* 
maker).  Secondly,  polynomials  can  never  be  completely  trusted. 
of  data,  even  though  they  are  from  the  same  application,  must  he  val  id-- 
(plotted  and  examined),  and  most  likely  adjusted  through  a ‘rial  and  r- 
process,  before  the  new  polynomial  fit  can  be  used. 


Verifica* ion  of  In’  rpolation  Sys’ems 


lU.  It  is  important  to  note  that,  derivat  ive  fau.1  ts  ar-  '-.si.l  v 
looked,  especially  when  the  user  restricts  his  evaluation  of  - inter: 
t.ion  procedure  to  testing  its  ability  to  reproduce  a physical  effect  v 
specified  error  bounds.  As  we  have  demonstrated,  one  can  have  a p--r f-  -■ 
satisfactory  model  of  the  observed  behavior,  while  that  same  model  is 
valid  in  other  important  physical  characteristics.  It  is,  therefore, 
highly  desirable  to  examine  all  characteristics  of  the  physical  system 
which  the  interpolating  model  is  expected  to  reproduce. 
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LARI  II:  CUBIC  .".PUNE  INTERPOLATION 


Tin?  aforementioned  reasons  provide  motivation  to  look  for  a 
: . t-i.y-t  -use  interpolating  technique  that,  will  give  smooth,  predict- 
: e:  -ivi  r ai.  i have  c-  ntinuous  first  derivatives.  U.  S.  Army  Engineer 

•■.•a;  Kxp-  riment  .’tat  ion  activity  in  several  engineering  areas 

2 ^ 

at.  r that  the  cubic  spline  " (hereinafter  referred  to  as  spline)  has 
ai  .-a:.*  a.-  .•  as  an  in*"rpolating  method.  Several  important  spline  char- 
• .s'  i s are  outlined  below: 

a.  The  splin"  passes  tlirough  th"  data  points. 

1_.  The  splin'1  is  a piece-wise  outic.  The  da1  a points  mark  the 
points  of  transition  from  one  set  of  coefficients  to  the 
next. 

c.  Th<  first  and  second  derivatives  of  the  spline  are  continu- 
ous at  all  points. 

d.  The  spline  curve  looks  similar  to  that  drawn  by  a french 
curve  or  a mechanical  spl ine  (more  on  mechanical  or  physical 
spline':;  in  paragraphs  20  and  33). 

e.  Ad.’ustjnent  of  any  data  point  affects  the  entire  curve,  but 
the  effect  is  predictable. 

f.  The  spline  is  uniquely  defined  by  the  X and  Y coordinates 

of  each  data  point  and  either  'he  first,  or  second  derivative 
at  each  data  point. 


Mathematical  Formula*  ion 

If.  The  mathematical  spline  fiinetion  is  a piece-wise  third  degree 
1 cubic)  polynomial  passing  through  all  data  points  and  having  continuous 
first  and  second  derivatives. 

17.  The  spline  can  be  viewed  as  a set  of  cubic  equations,  one  equa- 
’ion  for  each  interval  between  successive  da' a points.  The  coel l icients 
of  the  cubic  equal. ions  are  such  that,  at  any  data  point,,  the  equation  for 
the  left  interval  will  yield  the  same  values  for  the  first  and  second  de- 
rivatives, respectively,  as  will  the  equation  for  the  right,  interval. 

18.  Given  a set  of  N data  points  for  which  the  coordinates  (X^,Y^) 
and  second  derivatives  (M.)  are  known  for  every  point  (i  = 1»  2,  3...N^> 


where  i is  such  that  X.  < x < X 


The  first  and  second  derivative 


19-  The  method  of  evaluating  the  M.  coefficients  is  presented  in 
Appendix  A,  along  with  a derivation  of  the  above  expressions.  A compact 
fitting  and  interpolating  FORTRAN  program  is  included  in  Appendix  B. 
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Discussion  of  Spline  Characteristics 


20.  The  spline  properties  discussed  in  paragraph  15  illuminate  the 
advantages  and  disadvantages  of  splines.  Properties  a,  b,  c,  and  d combine 
to  produce  an  accurate  curve  having  continuous  first  and  second  derivatives 
This,  an  mentioned  earlier,  is  very  useful  in  engineering  analysis  problems 
Properties  d and  e provide  a physically  based  insight  for  adjusting  data  to 
obtain  a better  curve.  In  fact,  the  mathematical  cubic  spline  can  be  de- 
rived from  the  theory  associated  with  the  deflected  shape  of  a weightless 

elastic  beam  constrained  at  particular  points. 

2 

21.  Oreville  noted  that  many  physical  systems  are  correctly  modeled 
by  cubic  or  quadratic  equations.  This  fact  accounts,  to  a large  extent, 
for  the  popularity  of  cubic  splines.  It  also  enables  us  to  obtain  a 
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is  presented  in 
is.  A compact 
T-ndix  1-. 


r;.  'it  with  •«  .'mall  numb-r  of  ('<*'»  point.'  when  an  interval  of  the 
non  ■ 1 r ■ :ond,  r third  • t . ■ ■ ■/ ' r . 

:■  . . ' .-  nr  :.t  ep;  1 at  i ■ n 

• * * .ne  same  set  of  data  used  Li  fig.  1. 

• . smoothness  >f  the  curve  between  data 

" r*  ' .mpor*  ant  ir.  th-  ' ntinuous  leriva- 
. Thi:  is  a mar 1 improvement  over 
r ; v 1 !■  • 1 (.  t:  m del  shown  in  fig.  1 . 
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Fig.  2.  Spline  ap- 
proximation of  stress- 
stra. n relation  shown 
f ig . 1 


vei  ■ jmo- 

11  'rate:  in  fig.  3 whi  ’h  sh-  ws  a set 
hid  a fifl  i m- 

: pline  ha  been  fit.  rhi  xamj  . 

w.'  t:  influence  of  a "tad"  data  value  on  both 

' ■ • : 1 in*  and  polynomial.  Not'-  how  the  ,-plin< 

••  :.is  to  minimize  the  influence  of  the  t ad  point 
. Th"  polynomial  on  the  other  hand  is  -om- 
; I'-tely  upset,  by  the  baa  paint. 

2h.  This  rather  extreme  example  is  intro- 
duced not  a:*  a practical  consideration  but 
rather  to  illustrate  the  oscillating  character 
of  polynomials  that  is  so  troublesome  in  curve 
fitting..  Discussion  of  least  square  polynomials 
is  provided  in  the  applications  section,  para- 
graphs • . -J*‘  . 

Disadvantage,:  of  splin-. 

2r.  Properties  b and  f (paragraph  15)  are 
sometimes  interpreted  as  disadvantages.  The 
spline  is  not  defined  as  a single  expression 
over  the  entire  range  of  the  data  points.  There- 
fore, to  evaluate  S(x)  one  must  first  find  the  two  adjacent  data  points 
iunding  x and  then  apply  t.he  equations  of  paragraph  1^-  This  results 
in  longer  computer  time  than  that  required  in  polynomial  l Hs-  However, 
•arch  time  can  be  minimized  by  use  of  efficient  algorithms. 
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. • . rh<  i.umi  • r ;oefficien1  led  t lefini  • ■ plii 

times  the  number  of  data  punts)  Is  too  1 arg>  for  scum  engineer;-. 

iti on8 . rhis  is  due  1 the  flct  1 I I • mputer  codi  1 1 • ■ ■ 
a.-  the  finite  element  :u.  t finite  diff'-r'-n  ••  • le.  ; r<  ju.r>  n * -l,t.  . 

• ragi  ■ lve  physical  For  ins t • . • 

ysis  of  soil •structure  interact!  wi  problea  requii  that 

soil  zones  be  considered.  If  the  nonline ai  - trail 

oils  is  modeled  by  splines,  then  for  each  II:-  ■ • r» ■ i: * layer 
store  parameters  for  a different  split.-  fu:  • . . . nus , th*  :.um: 

storage  locations  required  by  splines  w.ll  b<  greatly  in  1 

problems.  The  same  is  true  for  finiti  lif ferenc 
problems.  The  engineer , quit-  naturall  , is  Inclined  t the 

storage  for  larger  physical  problems  ratt.-  r than  to  introduce  what 
feel  is  an  unnecessary  elegance  in  the  stre. v, -strain  inter;  olati  n r 
Summary 

27.  Although  WES  experience  with  splines  is  limited,  it  .a 
the  conclusion  that  the  advantages  of  the  s;lino  ( • ntinu  >us  i ri  . a* 
smooth  curve,  dependability,  and  physical  insight  far  utweigh  the  : 
lems  of  increased  computer  memory  and  somewhat  longer  computation  tim 
Further,  through  the  skill  of  the  numerical  analyst  and  programmer, 
can  exercise  considerable  influence  in  adapting  the  formulation  and  • 
to  minimize  run  time  or  storage,  a point  discussed  in  the  application, 
section. 


RATING  CURVE 


; a:-.  ::i;  aiti.i  -ati  • of  sflihes  to  engineering  problems 


T'  section  iescribec  the  applioati  n -'ll'  cubic  splines  to  sev- 
■ ns  WES.  Applications  of  the  cubic  spline 

• 1 ' • • • ri  !nt  ir.g  techniques  to  hydraulics  problems,  soil-structure 

■ tion  pr  « tat  . and  seeps  problems  ar<  presented. 


T.f  z U;‘ 


Hydrau  ics  In  : ler.s 


. r.g  -ur.n  s 

. : s applicat  i n in  the  field  of  hydraulic  engineering  deals 

•at  Ln(  irv<  f]  i routing  program  authored  by  Mr.  E.  A. 

■■  Mississippi  Valley  Division.  Mr.  Graves  and  Mr.  J.  B. 
'n.k,  Jr.,  >f  VIES  are  currently  involved  in  the  application  of  spline 
. • is  i this  problem. 

Rating  curves  are  used  to  describe,  at  stated  points  (stations) 

■ a river,  the  flow  (discharge)  characteristics  as  a function  of  the 

f the  water  in  the  river  (stage).  The  flood  routing  program  pro- 
iuces  very  satisfactory  results , but  a great  deal  of  difficulty  is  experi- 
enced in  btaining  polyn  mial  fits  to  the  rating  curves  that  are  required 
as  input  to  the  program.  The  f 11  wing  figures  and  discussion  illustrate 
'he  improvement  brought  by  the  spline  interpolation  method  to  the  rating 


•urve  model, 

•1 . Spline  fit  to  ?0  points. 
Fig.  1*  shows  20  data  points  of  a rat- 
ing curve  with  a spline  fit  through 
*-  se  points.  Except  for  the  large 
number  of  points,  the  spline  repr- - 
sentation  was  considered  satisfac- 
t ry  until  the  derivative  shown  in 
fig.  5 was  examined.  It  was  charac- 
terized by  an  objectionable  lack  of 
smoothness . 

T2.  Reducing  number  of  points 


to  improve  derivative.  In  order  to 


20-POINT 

SPLINE 


DISCHARGE 

Fig.  4.  Spline  fit.  of  20-point 
rating  curve 


RATING  CURVE! 


SPUNE 


DISCHARGE  T! 

Fig.  5-  Derivative  t r DO- point 

......  S( 

spline  lit 

P< 

the  physical  insight  to  the  spline 
has  its  impact. 

33-  As  mentioned  previously, 
the  mathematical  spline  formulation 
also  describes  the  deflection  of  a 
weightless  linear  elastic  beam  (a  phys- 
ical spline)  which  is  simply  con- 
strained at  the  data  points.  This 
means  that  we  can  visualize  the  mathe- 
matical spline  fit  as  that  shape  as- 
sumed by  a thin  strip  of  steel  which 
is  constrained  on  the  paper  by  a 
straight  pin  at  each  data  point.  With 


reduce  the  number  f data  [ 
pr  luce  a sr  ther  icriva’ ! < 

1 if  20  rigi.ial  p-  :.nt;  w*-r« 
and  a spline  was  f i + * > • d • 

This  is  illustrated  in  fig.  • 
with  the  original  ?0-f  In* 

. -point  urv  i ■ • • • 
riginal  to  th<-  desir'  d i<-.'r*-  . 
seems  desirable  to  add  a 
points.  In  a situati  n su 


RATING  CUR'/E 


POINT  w 

\< ' 

//  DO -POINT 


discharge 

Fig.  6.  Comparison  f • - a:,  i 
point  spline  fits 


20- POINT 


DISCHARGE 


this  concept  in  mini,  the  of:-  ■' 
adding  point  } (fig.  • ) to  the  - 
point  data  set  can  be  easily  pre- 
dicted. The  new  point  P would  t : 
the  spring  closer  to  the  original 
curve  along  A and  would  also  cau  ; 
a downward  movement  along  B. 

3^-  Fig.  7,  which  shows  the 
/-point  spline  fit  along  with  the 
original  20- point  fit,  shows  the  ef- 


Fig.  7.  Comparison  of  7-  and  20- 

point  spline  fits  Ct>  01  addlnS  point  P.  Note  how 


• ' ft)-  -urv>’.  By  some  fur 

ir. • ..1  i injjin  ved,  but,  it 

• th<-  j ii.*  ha:  been  made, 
t.).e  :-itt-  ■ Ti:  sj  line  makes 


r a i.ius'  '.ent 


ei;jen“  in  s::;  -othnes 


'Vement  i s 


;,:ei  nu:  • er  f pints.  This  asp*  - -t 
sj  ' ' interpolation  will  be  iis- 
assed  further  in  paragraphs  . 

- . J-’.  iil'i  :ati<  n of  spline  fi 


- ■musi-  rating  (or  conveyance  -urv 


I : . . I : . 1 . . : . I iVI  i ’•  i , . r r 

ujtaCH  ARGE. 

fc  ia  etimes  necessary  t l an<  - fig.  8.  Deri  itive  f r 7-point 

them.  With  splines,  this  can  be  done  spline  fit 

■y  subs1  itutii.,-  new  lata  p ints , with  the  assurance  that  the  new  fit  will 
e satisfactory.  This  coni  rai  • is<  f p lyn  rials,  which  re- 

ps i r*-  many  oore  data  points  and  which  must  be  carefully  studied  to  deter- 
" i ne  i:  a satisfactory  fit  has  been  btained.  Als-  •,  in  the-  application 
discussed,  the  first  derivative  is  r-  paired,  and  this  -an  be  obtained  di- 
r*-'tly  during  the  ep  ine  interp'  .1  at. i n pro  •>  dure.  Even  in  applications 
where  the  first  derivative  is  net  used,  if  a polynomial  curve  fit  is  to 
be  used,  it  is  desirable  t<  bfairi  tht-  first  derivatives  for  use  in  study- 
ing- the  suitability  of  the  --urve  t<  the  purp  in  hand. 
del  storage  curves 

37-  The  success  with  modeling  rat  ing  curves  soon  led  to  experimen- 
tation with  modeling-  storage  curves.  This  flood  r<  utinp  application  has 
severe  computer  menr>ry  restrictions  n it,  so  extra  steps  were  taken  to 
reduce  the  coefficients  stored.  This  was  accomplished  by  retaining  only 
the  coordinates  of  the  points  and  beam  moment  of  each  data  point.  his  is 


original 
ilso  ause 


shows  the 


with  the 
lows  the  ef 


Note  h w 


restricted  memory  envi  r iirr.ent,  j nving  1 r t.h*'  modi  :'i  <-nti  r witli  • -r* 
l nterpo latino  coraputati  -n  ' in’  . 

Phe  f rn ilatl  i ale  pr  /idei  inear  extra]  Latj  n 1 ■ 
bey  nd  the  lata,  us  Lng  th<  p<  and  • rdj  Ltes  i h<  earesl 

] int . (A  tabulati  n f this  pr  ra  pi  ented  Li  Appendij 

r:  dification  gives  an  ext  ra  advantage  w>  • n used  in  a • nvergin,'  • r : 
error  procedure,  the  advantage  h«-i.  * ;.v  * • . • ria..- 

side  the  rang*  of  the  physical  data.  ihe  sf  i no  • rrruia*  : mak*  ' 

to  assure  that  extrapolated  values  wi  aid  ■<  :..  < rg<-ri  ••  : • .<•  • ria 

error  process.  Such  assurance  is  btain<  1 oi  I iditioru 

ming  in  polynomial  fits.  It  is  a sufficiently  Jiff 
polynomial  within  the  bounds  of  the  lata  wit:,  ut  • • : < ring  with  it 
the  physical  data. 


1 i : -St.ru  :tur<  Inters  ■ 1 i >n  Studies 


39-  In  soil-structure  inters  :ti<  n problems  using  the  finite  differ- 
ence or  finite  element  computer  codes,  it  is  necessary  to  represent  t:.< 
nonlinear  stress-strain  behavior  of  the  s Us.  is  i.;  normally  inpu*  as 
sets  of  coefficients  for  polynomials , one  set  for  each  soil.  Using  this 
approach,  it  is  possible  to  compute  the  shear  and  bulk  moduli  of  the  ma- 
terial at  any  stress  or  strain  level . The  typical  solution  procedure  re- 
quires the  computation  ol  the  first  derivative  of  the  curve  in  addition 
to  evaluating  functional  values  from  the  -urve  itself.  Spline  fits  are 
ideally  suited  to  represent  the  curves  because  the  fit  is  smooth,  the  de- 
rivatives are  continuous,  and  only  an  intermediate  skill  level  is  needed 
to  obtain  a satisfactory  fit  to  the  data. 

Spline  versus  polynomial  fit 

40.  The  89  data  points  shown  in  fig.  9 are  an  oxample  of  those 
modeled  in  programs  of  this  type.  The  points  were  computed  directly  from 
measurements  taken  during  a test  of  a soil  specimen.  Both  a spline  fit 
and  fifth  degree  least-squares  polynomial  fit  are  shown.  The  two  points 
to  the  left  of  the  pressure  axis  were  added  to  force  the  polynomial  to 
pass  through  the  (x=0,y=0)  point. 
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•■1 . Judgi  ng  n ■ y f roi  fig.  • 

uld  easily  conclude  that,  the 


5 DEGREE 
POLYNOMIAL  . 


SPLINE 


FIRST  INVARIANT  (X) 


Reducing  number  of 
)oints  to  improve  derivative 

42.  Does  the  physical  effect 
being  modeled  require  all  89  points 
for  its  correct  definition?  Inspec- 
tion of  fig.  9 indicates  this  is  un- 
likely. The  effect  on  the  spline  of 
the  errors  in  the  data,  small  though 
they  may  be,  must  be  considered. 

Note  that  the  effect  of  any  data 
error  is  magnified  as  the  points 
are  brought  closer  together . The 
large  changes  in  slope  shown  in 
fig.  10  indicate  that  the  spline 
is  being  subjected  to  many  torques 

and  countertorques  that  are  principally  due  to  the  data  error  and  the  close 
spacing  of  the  points. 

43.  There  are  smoothing  spline  methods  that  overcome  this  problem 
(by  not  passing  through  all  data  points).  Unfortunately  they  require  re- 
peated trials  and  evaluations  to  obtain  a satislactory  value  for  the 


SPLINE 


5 DEGREE  POLYNOMIAL 


additi 


FIRST  INVARIANT  (X) 


Fig.  10.  Comparison  of  spline  and 
polynomial  derivatives  for  89-point 
fits 


smoothing  parameter.  From  a practical  star.  Ij  int,  it.  in  geri'T's  1 . pr- 
able  and  less  time  consuming  to  select  a few  lata  points  a*  ').•  u* . < • 

rather  than  to  make  repeated  runs  in  search  f a g i parai-  •'  er  r 

new  set  of  data.  (This  is  not  to  imply  that  ne  method  is  sup'  r:  r ' 
the  otinr,  rather  that  interpolating  splitv-r  can  <•—  r.  mi  a .y  h<  mad-  • 
function  satisfactorily  for  many  scientific-engineering  ajj^i-ati  ns. 

44.  Fig.  11  shows  spline  and  polynomial  fits  1 L<  f th<  ri 

89  points  of  fig.  o.  N te  t • a* 
] <ints  shewn  are  not  the  resu..  * 


SPLINE 


several  trials  to  obtain  the  t 
set . 00  , wer<  singly  selecl 

that  the  x values  were  about  n 


00  5™  DEGREE  * inch  apart  on  the  original  dra* ' 

uJ  [ /POLYNOMIAL  ^ . -•*"  a procedure  which  certainly  car.: 

^ be  called  "tuning"  the  data. 

45.  Fig.  12  shows  the  vast 

1 — > — > — > — 1 — 1 — * — 1 — 1 1 — * j 1 

FIRST  INVARIANT  /XI  1",prove'i  spUne  ltri”“i"  plot 

which  now  appears  superior  t 1 ' 

Fig.  11.  Comparisons  of  lf-p  Int  the  89-  and  16-point  polynomial 
spline  and  polynomial  fits 

lenvatives  . Any  improvements  r<  - 

quired  can  be  obtained  with  the  aid  f the  aforementioned  insight  about 
the  physical  equivalent  of  the  cubic 


POLYNOMIAL 


FIRST  INVARIANT  (X) 

Fig.  11.  Comparisons  of  l6-]»  int 
spline  and  polynomial  fits 


spline  (paragraph  33)- 
Summary 

46.  The  authors  have  been  di- 
rectly involved  in  three  soil- 
structure  interaction  projects  at 


SPLINE 


DEGREE  POLYNOMIAL 


FIRST  INVARIANT  (X) 


R.  E.  Walker  and  fig-  12.  Comparison  of  spline  and 

Cheek  (March  1970)  on  polynomial  derivatives  for  16- 
a dynamic  nonlinear  point  fits 

elastic  finite  element  method  stress  analysis  program. 

J.  L.  Kirkland  and  Cheek  (April  1970)  on  a nonlinear  elastic 
incremental  loading  soil-structure  interaction  program. 

PFC  B.  Phillips  and  Radhakrishnan  (1970)  on  a finite  differ- 
ence solution  of  nonlinear  elastic  analysis  of  ground  motion. 
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•«  ; :i-  . 1 . ■ V • 1 V !.'•  -I  1 u ; i ' J I . ; I . • ' uo  ..  1- 

•u  • ;r<  it.’  • rn  •*  i,  • n . is  w rthwhili  . However,  it  ’an not  be  stated 
rtainl  1 t thi  puted  results  r better,  altl  igh  there  an 
...  ement  in  tl  stabilit;  f tl  proc<  iur-  ( :■  • - redue- 

* ra.t.  energy  growth),  since  no  known  correct  solution  with  which 
• result!  cists.  It  is  felt  that  splines  have  removed  some 
rrors  ii  l<  nlinear  characteristics  and  1 ave,  at  the 

• . fed  us  t 'a-  tr  r<  sear  :h  ther  probl  ems  • 


. . Producing  a contour  map  from  a 


atisfact  rily  acc  i ; 1 i shed . A special  case  of 
occur  on  a rectangular  grid  has,  how- 
Kig.  13  shows  the  data  configuration. 

GRID 

RECTANGLE 


has  not  in  general  been 
':.is  probl e.-:  where  the  data  p hit 
■cor,  been  . ved  us  i tig  splines. 


COL  MN  . 


at  ur  map,  a surface 
■ ust  first  in  fitted  to 


this  seethe,  a methc  i ... 

, ♦ ♦ ♦ 

which  fits  ■!  simple 

, * — 1 — 

bicubic  spline'*  to  the  •'lO'ti  ' ^ ! 

data  is  described,  and  ... 

an  application  to  con- 

fined  steady- state  ..  - 

Fig.  13-  Data  configuration  for  spline 

seepiage  under  a weir  contouring 

is  discussed. 

ho.  The  first  part  of  the  procedure  is  to  1 it  a cubic  spline  al 
<ach  row  and  each  column  of  the  data  points 


* 

D~/ 

• ♦ 

A 4 

♦ • 

• * 

♦ ♦ ♦ 

V 1 
4 * 4 

4 H 

♦ ♦ 

1 ‘ 

i * * 




^here 


Z.  = Z value  c rr^-s ponding  to 
1 >J 

The  cubic  spline  as  defined  n the  j 
property  f f paragraph  15  is 


l,j  is  the  second  derivative  of  the  vertical 


through 


Krom  these  cubic  splines,  the  simple  bicubic  spline  is  generated 
on  (grid  rectangle) 


,_v  is  t he  bicubic  spline  as  defined  in  the  grid  rectangle/ 
ft-hand  c rdinates  are  (X.,Y^)  . From  this  bicubic  spline 
instruction  of  the  contour  lines  can  be  accomplished  as  ex- 


it^ PERVIOUS 
' STRUCTURE 


MPERVIOUS 


Fig.  14.  Steady- state  seepage  under  a weir 


ti  c distribution  >f  pressures  in  the  pervious  medium.  Th<  result 
n rmally  expressed  as  contours  of  equai  ] ressures  r potentials 
equipotentials . 

52.  Region  ABCF  in  fig.  14  is  divided  into  a grid  sys 
finite  difference  solution  to  tiie  governing  partial  differentia: 
(Laplace's  equation) 


ith  boundary  conditions 


0 on  ED  and  AB  (ED  and  AD  are  flow  line, 
flow  across  these  boundaries.) 


There  is 


the  output  I rom  this  solution  is  used  as  input  to  the  - 

>grair  to  obtain  the  equi potential  lines . The  results  are  sh 
The  equipotentials  are  labelled  as  the  percentage  £ hp  - h, )/ 
100  , where  h(  is  a given  equipotential . 

.he  quality  of  the  contour  map  can  be  illustrated  in  three  ways 

a.  The  contour  lines  are  smooth. 

b.  The  problem  was  set  up  such  that  the  center  of  the  weir  was 

r,f  the  region.  One  would  therefore  expect  th 
equipotential  line  to  be  a line  of  symmetry.  This  is 
indeed  the  case. 

c.  The  boundary  condition  dh/dy  =0  on  ED  and  AB  require 
that  the  equipotential  lines  intersect  the  line  segments 
orthogonally.  Again,  this  is,  in  fact,  the  situation. 


U . It  is  concluded  that  the  simple  bicubic  spline  is  an  excellent 
function  for  contouring  data  established  on  a grid. 
llnconfined  seepage 

55-  A problem  of  classic  importance  in  civil  engineering  analysis 
is  that  of  unconfined  flow  in  an  earth  darn,  with  special  interest  given  to 
the  phreatic  surface, ^ ' which  is  a top  flow  line  across  which  there  is  no 
flow.  Another  property  of  the  phreatic  surface  is  that  the  potential  at 


PHREATIC 

SURFACE 


FINITE  DIFFERENCE 
GRID 
7 H W 


SURFACE  OF 
SEEPAGE 


Location  of  phreatic  surface 


with  boundary  ■■  ndit.ion, 
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‘he  procedures  used  in  s<- Lying  this  pr-biem  are  as  follows: 

2-  Make  an  initial  gue.n;  to  the  phreatic  surface  >:  . 

j_  • Solve  the  pr  h'em  by  finil*  iifferonce  retd!  is  as  . : 
a confined  seepage  problem  using  db/dn  -ml  ; /- 

stani  ' pi  isur  : p ■ ass  iensil  • i 

g ai-celerat ion  duo  t gravity). 

£.  Adjust  DC  t.o  satisfy  h y . 

d.  Repeat.  £ and  _ until  dh/dn  = 0 anl  h r,r«  . -r 

on  DC  simultaneously. 

! • The  principa I llfficult  f thl  pr  :edure  is  in  s&tisl  : 

Sh/dn  = 0 on  DC  . Fhls  is  becausi  intersects  the  finiti  • ' • • 

: . i .it.  irregular  t Lnts,  as  is  t rat  < 1 . . . ! r.  ■ •••  ratii 

• th<  I r ■■  lure  may  e ; a*. i this  LiffJ  . ty.  ; l 

should  then  be  altered  as  foil  ws : 

a.  Fit  a cubic  8|  in<  al  >ng  DC  as  shown  in  fig.  lb.  Note 
that,  the  points  of  intersection  of  the  phreatic  surface 
and  the  grid  are  used  as  data  points.  Use  the  set  of  equa- 
tions which  solves  for  the  slopes  at  the  data  points  (see 
Appendix  A for  details),  and  set 

S = -cot  a 


SM  = -tan  0 

where  N is  the  number  of  data  points. 

Replace  3h/&n  = 0 and  l/pg  = constant  on  DC  with  the 
equivalent  expressions 


m,  - 


1 + 52 

l 


(*1)  - 3i 

1 , G‘° 


1,  2,  3...N 


where  the  S^s  are  the  slopes  as  computed  by  the  spline 

I lh  • 

c.  Incorporate  these  simple  formuLas  into  the  finite  difference 
solution  to  obtain  the  next  trial  values. 


F/.i-.T  IV:  CONCWr.JONf 


The  spline  interpolation  technique  is  a valuable  tool  in  curve 
fitting  l’or  oumputer  programs  for  which  more  commonly  used  techniques  are 
unsuitable  or  of  limit  i value.  The  spline  technique  has  three  primary 
advantages : 

a.  The  spline  passes  through  all  data  points. 

b.  The  first  and  second  derivatives  of  the  spline  are  continuous 
at  all  points. 

£.  The  spline  can  be  easily  modified  to  satisfy  new  or  addi- 
tional data. 

WES  experience  in  applying  spline  techniques  to  engineering  problems  has 
indicated  that  in  many  applications  these  advantages  outweigh  the  addi- 
tional storage  and/or  computation  time  requirements  of  the  technique. 

60 . Since  the  spline  function  is  required  to  pass  through  all  data 
points,  erratic  derivative  behavior  may  result  from  experimental  error 
when  the  data  points  are  numerous  and  closely  spaced.  Trial  and  error 
methods  for  smoothing  such  functions  exist,  but  they  are  time  consuming. 

WES  experience  has  indicated  that  acceptable  results  can  generally  be  ob 
tained  by  simply  selecting  a more  limited,  more  widely  spaced  set  of  the 
data  points  to  which  to  fit  the  curve. 
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APPENT1X  A:  EQUATION:;  FOH  CUBIC  CPLlIfE 

1.  A cubic  spline  is  a function  whose  third  derivative  is  a step 
fun-tion  with  points  of  discontinuity  a*  the  data  points 


Xi’Yi 


i = 1,  2,  3...N 


where  N is  the  number  of  data  points.  The  spline  function  of  degree 
three  is  therefore  a piece-wise,  continuous  third  degree  polynomial  having 
continuous  first  and  second  derivatives. 

2.  For  this  application  the  function  will  be  further  required  to 
interpolate  the  data  points.  Hence,  on  the  i interval 


X.  i x s 
l l+l 


i = 1,  2,  3...H  - 1 


S(x)  = a.  + b.x  + c.x2  + d.  x3  (A1 

ii  i i 

where  a , b.  , c.  , and  d.  are  constants  to  be  evaluated.  Equation 

i’ll  i 

A1  can  be  rewritten  as 

S(x)  = a/(x  - Xi)3  + ^(Xi+]_  - x)3 

+ c^(x  - x±)  + d^(xi+1  - x)  (AS 

where  a'  , b'  cJ  and  d/  are  an  alternate  set  of  constants.  Dit- 

i ’ i i i 

ferentiating  equation  A2  twice  yields 

S"(x)  = 6a/ (x  - X±)  + 6b'(Xi+1  - x)  (A/ 


Applying 


S"(X. ) = M. 

v l l 


S"<W  - Mitl 


AI 


where  Nh  is  the  second  derivative  at  point  (x  ,Y, ) , the  result  is 


and  b ' — . i , 

i MX.  - X, 
i +-3  i 


3*  Since  the  splice  nual  interpolate  tin  lata  pointi 


< m>  rj 


ubstitutinc  these  into  equation  A2  produces 


M / 

y . (x  . , 

li+i  Z \ ‘i+i  1 


X - X 
i+1  "i 


Collecting  terms 


- 1 ..  (xm  ~ *)3  1 (*-*i)3 

hn-h  s ln  Vr7! 


Xi+1  " x 
xw  - Xi 


Y v \ 2 (X  Xi) 

Xi+1  xi ) x TT 

/ J Ai+1  Ai 


(A8) 


/ x • , , - X 

S"(x)  - M.  — — I + M. 


A . 
1 


i l X,xl  - X. 
V i+3  i, 


'"'Vi  - 


(A13) 


xi  * x * Xi+l 


i - 1,  2,  3 - - -n  - 1 


(A9) 


U.  The  NL's  are  evaluated  by  satisfying  continuity  of  the  first 
rivative  at  the  data  points 


X.  , Y . 

l l 


2,  ■ 


•uid  speeifyin  • M,  and  M . At  point  i 


(AID) 


S'(X.-)  --  : ’(X.  + ) 


A n lying  equation  A12 

M.  Y,  - Y . - 

(x,  - X,  ,)  + -/  . *■  7 (M,. 


t - w - 1 ■■  ,-:y  ♦ 7 («i-i  - Hi)(xi  - t-d 


(All) 


- T <xi,i  - xi>  * V,',;-.  ' z (Mi  * Mi+i)<xi*i  * xi>  (JU*> 


Collecting  terms 


A3 


69647 


A good  choice  for  M,  and  M is 

1 H 


This  set  of  simultaneous  linear  equations  can  be  solved  for  the  M's  . 

5*  H is  sometimes  more  desirable  to  solve  for  first  derivatives 
at  the  data  points  rather  than  second  derivatives.  This  set  of  simul*a- 
neous  equations  is  derived  as  follows.  Satis tying 


where  T^  is  the  first  derivative  at  (X.,Yj ),  equation  A12  becomes 

Ti  ' - T (xin  - V ♦ 777  . x'  * Z <Mi  - Min,(xin  ' xi>  <u?> 


«i  ■ “ 0 


(AI7) 


(Alb) 


Y,.,  - Y. 


W - -Y  <xin  - V ♦ - 


Collecting  terms 


Ci  ' (■  ^ - V)  (Xln  - X1 
'in  * * “T1)  (xm  - xi 


Yj+1  " Yi 

Xi+1  ‘ Xi 


(xm  - V * T~rr 

i+1  i 


(Y,.,  - Y«  Y j 

-i  + _ L 

X...  - x.  x. 

i+l  i i 


(A21) 


(A22) 


Solving  for  M.  and  M.  , , 
i l+l 


(A15) 


Mi  = 3X7 


Y.  . . - Y. 
l+l  i 


i+1  ' Xi 


' Ti+1  “ 


...  - X. 
l+l  l, 


(A23) 


v _ y \ 

M = IT.  + 2T  - 3 )f- — v 

!+i  v 1 1+1  x.+1  - x.Axi+1  - Xiy 


(A2U) 


(Alb) 


'i  s . 
•fives 
multa- 


(A17) 

(A16) 


Requiring  continuity  of  the  second  derivative  at  data  points  i = 2,  3> 

- 1 


SH(Xi+)  = s"(x.-) 


(A25) 


Applying  equations  A12  and  A2l 


Yi+1  ’ 


Xi+1  ' X 


= T.  + 2T. 
i-l  l 


A - yi-iV__£_ 

xi  - Xi-Jlxi  - xi-i 


(A26) 


Or 


(xi.i  - xi)Ti-i  * 2<xi+i  - xi-i)Ti  * <xi  - xi-i)Tiu 


) (A19) 


±)  (A20) 


= 3 


X_,  - X.\  /x  - X 

11  : (*i  - w * <Vi  - V 


ixi  - xi_i, 


i+i  iy 


(A27) 


i = 2,  3,  4...N  - 1 


(A21) 


(A22) 


M = = 0 becomes 

1 N 


Y - Y 
2 1 

2T  + T - ^ — = 

* 1 2 J x2  ~ xx 


T + 2T  = 3 
N-l  N J 


yn  - yn-i 

XN  ' XN-1 


(A28) 

(A29) 


This  system  may  be  solved  for  the  first  derivatives  at  the  data  points, 


At) 


IX  B:  : re  : n ING  AM  INTER!  LAI  ENG  SU1  UTINES 


his  appendix  ‘ont&ine  listings  f the  two  FORTRAN  language  sub* 
stines,  SI  INI  (table  B1‘  ) and  SPLINT  (table  B2).  The  following  sections 
. 'rib'-  th*  .se  of  the  two  routines . 

, ubroutine  i LINE 

?.  Th'  . ! LINE  subroutine  is  used  to  fit  a cubic  spline  to  a set  of 
ix,y)  data.  That  is,  it  calculates  the  moments  at  each  ii  terior  data 
pc int,  assigning  zero  to  the  moments  at  the  end  points.  Note  that  only 
three  arrays  are  required  for  each  spline,  an  x , ay,  and  a moment 
array.  The  spline  fit  is  accomplished  by  a call  statement  in  the  user  s 
program  of  the  form: 

CALL  SPLINE  (Al,  A2,  N3,  Ah) 


where  the  arguments  Al,  A 2,  N3,  and  Ah  are  as  follows: 


'•.rpumcnt 

Al 

A2 

N3 

Ah 


Purpose  

Names  the  independent  variable  array  (x). 

Names  the  dependent  variable  array  (y). 

Specifies  the  number  of  (x,y)  points.  N3  must  be 
in  integer  form. 

Names  the  array  into  which  SPLINE  will  store  the 
calculated  moments . 


3.  Note  that  the  user  must  specify  the  size  of  arrays  Al,  A2,  and 
Ah  through  DIMENSION  or  COMMON  statements.  Several  splines  may  be  fit 
and  saved  for  subsequent  use  by  making  successive  calls  to  SPLINE  using 
different  names  for  arguments  Al,  A2,  and  Ah  and  indicating  the  number 
of  points  through  argument  N3. 

h.  Example:  Given  two  sets  of  (x,y)  data,  fit  a spline  to  each 

set  of  data. 


B1 


‘ . Solution:  Let  one  set  of  data  be  in  arrays  XI  and  Y1  having 
Nl  points.  Let  the  other  set  of  data  be  in  X2,  Y2  having  N2  points, 
fol lowing  FORTRAN  statements  are  required,  assuming  that  there  are  no  : 
than  30  points  in  either  lata  set. 

DIMENSION  Xl( 30),  Yl(30),  Cl(30),  X2(30),  Y2(30),  02(3°) 

Statements  to  input  the  two  data 
sets  and  specity  their  size  in 
NI  and  N2.  (Note  that  Cl  and  C2 
need  not  be  set  to  any  specific 
value . ) 

CALL  SPLINE  (XI,  Yl,  Nl,  Cl) 

CALL  SPLINE  (X2,  Y2,  N2,  C2) 

' . At  this  point  array'  Cl  will  contain  Nl  moments,  and  array  C2 
will  contain  N2  moments;  the  first  and  last  moment  in  each  array  will  be 
zero. 

7 . Note  that  Nl  need  not  equal  N2  but  neither  may  be  less  than  2 
nor  greater  than  the  maximum  size  specified  for  the  associated  arrays  ( 3 
in  this  example). 

Cubrout ine  SPLINT 

L.  Subroutine  : i LINT  (SPLine  INTerpolate)  operates  on  a cubic  spline 
delined  by  the  (x,y)  coordinates  of  N points  and  the  moment  at  each  of 
those  points.  It  calculates  values  of  y and  y’  for  any  value  XX  of 
the  independent  variable  x . Should  the  value  XX  lie  beyond  the  range 
01  data  defined  by  the  points,  this  program  will  extrapolate  linearly  from 
the  first  (or  last)  using  the  slope  of  the  spline  at  that  end  point. 

9.  Interpolating  is  accomplished  by  a call  statement  of  the  form: 

CALL  SPLINT  (Al,  A2,  A3,  AL,  A5,  A6,  N7) 
where  the  arguments  A1-A6  and  N7  are  as  follows: 


irgument 


Purpose 

A1  Specifies  the  value  XX  of  the  independent  variable 

x for  which  y and  y'  are  desired. 

A2  Receives  the  vtilue  of  y at  XX  computed  by 

SPLINT. 

A ^ Receives  the  value  of  y'  at  XX  computed  by 

SPLINT . 

A4  Names  the  independent  variable  array  of  the  spline. 

Ar:  Names  the  dependent  variable  array  of  the  spline. 

A6  Names  the  moment  array  of  the  spline. 

N7  Specifies  the  number  of  (x,y)  points  that  define 

the  spline.  N7  must  be  in  integer  form. 

10.  Note  that  arrays  At,  At,  and  A6  must  contain  N7  values  each  of 
x , y , and  moment,  respectively;  i.e.,  they  contain  the  spline  defining 
data.  Argument  A1  is  an  input  argument  for  XX  , while  A2  and  A3  receive 
the  values  for  y and  y ’ calculated  by  SPLINT . 

11.  Two  of  tue  variables  in  subroutine  SPLINT  may  be  useful  in  some 
applications.  Variable  FPPXX  contains  the  value  for  y"  . Variable  M 
indicates  whether  the  computation  was  an  interpolation,  in  which  case  M 
is  zero,  or  an  extrapolation,  in  which  case  M is  -1.  A reduction  in  run 
time  would  likely  result  from  incorporating  a more  sophisticated  search  pro- 

edure  than  that  used  ( statement  numbers  100  through  140  in  table  B2) . 

12.  Example : Assuming,  that  the  steps  outlined  in  the  example  for 

subroutine  SPLINE  have  been  taken,  the  following  statements  would  calculate 
y values  (YY)  and  y’  (YPRIME)  at  XX  = 36.49  from  the  second  set  of 
spline  data. 

xx  = 36.49 

CALL  SPLINT  (XX,  YY,  YPRIME,  X2 , Y2,  C2,  N2) 

Test  Program 

13.  Table  B3  shows  a simple  test  program  and  the  ten  interpolated 
values  computed  using  the  GE  430  Time  Sharing  system  at  the  U.  S.  Army  En- 
gineer Waterways  Experiment  Station. 

B3 


Table  B1 

Subroutine  SHJ’JE 


1 000 

SUBROUTINE  SPLINE  (X,  ZY,  N»  52) 

1 Oioc 

SPLINE  FI  TUNG  SUBROUTINE  ADAPTED 

f-nJM 

1 020C 

WORK  BY  GREVILLE,  U 5 AhMY  MA1H. 

HESEAhCH 

1 030C 

UN  IV  Of-  WI5C0NS0N,  T 5 REPORT  \893. 

1 040C 

JUNE  1968. 

lObO 

DIMENSION  X(l),  ZY(P»  52(1) 

1 060 

DATA  EPSLN  /l.L-6/ 

1 070 

N1  = N - 1 

1080 

ASSIGN  110  TO  I SW 

10  'XI 

DO  130  1 = 1 . N 1 

1 100 

H = X(I  + n - X ( I ) 

1 1 10 

DLY  = ( Z Y ( I + 1)  - ZY( I ) ) / H 

1 IPO 

GO  TO  I SW,  C110,  100) 

1 130 

100 

HPZZZ  = HL  + H 

1 1/10 

52(1)  = 2.  * C DLY  - YL)  / H2ZZZ 

1 150 

GO  TO  120 

1 160 

1 10 

ASSIGN  100  10  ISW 

11  70 

120 

HL  = H 

1 180 

YL  = DLY 

1 1 90 

1 30 

CONTINUE 

1 POOC 

1 P10 

S2( 1 ) = 0. 

1 PPO 

S2CN)  = 0. 

1P30 

OMEGA  = - 1 .0717968 

1 240 

140 

ETA  = 0* 

1 250 

ASSIGN  170  TO  ISW1 

1 260 

DO  190  I * 1,  N 1 

1 2 70 

H * XCI  «•  1)  - XCI) 

1 280 

DLY  = CZYC I ♦ 1 ) - ZYC I ) ) / H 

1 290 

GO  TO  ISW1,  (170,  150) 

1 300 

150 

H2ZZZ  = HL  ♦ H 

1 310 

HI  = .5  ♦ HL  / HPZZZ 

1 320 

W = (HI  ♦ S2 ( 1 -1)+  (.5  - BI ) * 

S2( I ♦ 1 

1330% 

3.  * (YL  - DLY)  / H2ZZZ ) * OMEGA 

1 340 

S2(  I ) = S2(  I ) + W 

1 350 

Z = ABS(W) 

1 360 

IF  (Z  - ETA)  180,  180,  160 

1370 

160 

ETA  = 7. 

1 380 

BETA  = S2C I ) - W 

1 390 

GO  TO  180 

1 400 

170 

ASSIGN  150  TO  ISW1 

1 410 

180 

HL  = H 

1 420 

YL  = DLY 

1 430 

1 90 

CONTINUE 

1 440 

IF  ( ABS( BETA ) * EPSLN  - ETA)  140, 

140,  200 

1 450 

200 

CONTINUE 

1 460 

HE TURN 

1470 

END 

GEN  1 


52(  I ) ♦ 


Table  b 2 

Subroutine  .T! 


1 520 

SUBHOUIINF  SPLINi  (XB.  FXX,  FPXX,  X, 

ZY  , 

S2.  N) 

1530C 

3PLINF  INIFkPOLAI I NG  SUBHOUIINF  » bY 

DAY 

C H F t K 

1 540 

DI  MFNSI  ON  X(l),  ZY(1).  32(1) 

1 550 

MM  = 0 

1 560 

XP  = XL) 

l 570 

1 = 1 

l 580 

IF  (XP  - X(  1 ) ) 100/  170, 

1 10 

1 590 

1 00 

MM  = - 1 

1 600 

XP  = X(  1 ) 

1 610 

00  10  170 

1 620 

1 10 

IF  ( XP  X ( N ) ) 130,  150, 

1 40 

: 630 

1 20 

lb  ( XP  - X( 1 ) ) 160,  170. 

1 30 

1640 

1 30 

I = I 1 

1 650 

GO  10  120 

1 660 

1 40 

MM  = - 1 

1 670 

XP  = X(N> 

1 68  0 

1 50 

1 = N 

1 690 

1 60 

1 = 1-1 

1 700 

170 

Hi  1 = XP  - X(  I) 

1710 

H12  = XP  - X( 1 ♦ 1 > 

1 720 

PhOL)  = HU  * Hi  2 

1 730 

L)X  = X(I  ♦ M - X ( I > 

1 7 43 

DELY  = (Z  Y(  I ♦ 1 ) - ZY<  I ) 

) / DX 

1 750 

S3  = ( S2(  1 ♦ 1 > - S2( I > > 

/ DX 

1760 
1 770 
17B0 
1 790 

FPPXX  = S2( I ' ♦ HI  1 * S3 
UELSQS  = (S2(I)  ♦ S2( I + 
FXX  = ZY(  l ) ♦ Hit  * DELY 
FPXX  = DELY  + (HU  + HI  2) 

1 ) ♦ FPPXX)  / 6. 
•f  PROP  * DFLSQS 
* DELSOS  ♦ PhOD 

* S3  / 6 

1R00 

IF  (MM.EO.O)  GO  TO  180 

18  10 

FXX  = FXX  ♦ FPXX  ♦ (Xb  - 

XP) 

1 820 

1 80 

CONI INUF 

1 830 

KETUKN 

1 840 

END 

HF  A DY 


I 


Table  B3 

Test  } r .gram  f’  r .TT  ! ar.3  ■’!  1 I '■”!  ■'  -D  r if  ■ ^ _ 

10in  DIMENSION  X(im,  Y ( 10)  i C ( 10  5 

1OP0C  sKT  1 HI-  TES  1 DATA. 

1030  X<  1 ) = 1 -ft 

1 0 AO  Y<  1 ) = 1 . 

1050  X<2)  = 5. A 

1060  Y(2>  = 2. 

1070  X ( 3 ) = 7. 

1080  Y(3)  = 1. 

1090  X ( A ) = «.2 

1100  Y( A)  = 1 . 

1110  NUMB  = A 

1 120C 

1130C  FIT  A SPLINE  THhOUGH  THE  X,Y  LATA  POINiS. 

1 1 AO  CALL  SPLINE  (X,  Y » NUMB,  C) 

1 1 50C 

1160C  I NTEHPOL.ATE  AT  INTERVALS  OF  ONE. 

1 1 70  PHINT  300 

1 IPO  DO  100  1 = 1 , 10 

1190  X I = I 

1200  CALL  SPLINT  (XI,  YY,  YP,  X,  Y,  C,  NUMB) 

1210  100  PHINT  200,  XI,  YY.  YP 
1220  STOP 

1230  200  FORMAT  ( 3X  3E20.9) 

1 2A0  300  FORMAT  (//10X  1HX,  1 9X  1HY.  19X  7HY  PRIME  /> 
1250  END 

READY 

RUN 


OP  8 3 1 WES  05/19/71 


0. 1 OOOOOOOOE+O 1 
0 • 200000000E+0 1 
0 • 300000000E+0 1 
0. AOOOOOOOOE+Ol 
0 . 500000000E+0 1 
0 • 600000000E+0  1 
0 « 700000000E+0  1 
0 . 800000000E+0  1 
0 • 900000000 E*0 1 
0. 1 0000 0000 E *02 


RUNNING  TIMEt 


0. 605953327 E +00 
0. 126029A07E+01 
0. 1RAP63327E+01 
0.219698582E+01 
0.216050A1  3E-M11 
0. 16091 7 1 68  E ♦ 0 1 
0. 1 OOOOOOOOE  + O 1 
0.96708P5A6E+00 
0. 1 1 35431 8 1 E+0 1 
0. 1 30A721 5BE  + 01 


3.2  SECS  I/O  time 


Y PRIME 

0.6550777881 *00 
0.6A20A9980E+00 
0. A95A87 1 39h *00 
0. 186076697E*  0 

- .28618 1 3A6E*00 
- .727 1 31 570E+00 
- . 3385795301 *00 
0. I 551822P5E*  0 
O.  1 69289765E  + 00 
0. 1 692897651 *00 
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Til  IPONlONINO  MILITARY  ACTIVITY 


This  paper  wv  p>i"par<  i • familiarize  practicing  scientists  and  engineers  with  the 
ill.  1 ■ ; 1 inf  interpola’  : • te  -hnlquc  as  a possible  tool  in  curve  fitting  for  computer 
program:.  ! r wt  i ••>.  ».or>-  o-  l.ly  use  i techniques  may  be  unsuitable  or  of  limited 

it  i ■ re  ethods,  spec  i f ical  ly  piece- 

wire  linear  and  j lynomial,  and  examples  f applications  of  the  technique  to  engi- 
neering pi  blems  are  presented.  Appendix  A contains  the  mathematical  derivation  of 
’h«-  "tiiati'i.  defining  th.  spline  function,  and  Appendix  B contains  a compact  FOHTRAf; 
fitting  and  . n*  erp'  dating,  j r ,-ras..  The  interpolating  spline  curve-fitting  technique 
has  three  primary  advantages.  ( a)  the  spline  passes  through  all  lal  a points;  (b)  the 
first  and  ■ mi  derivatives  f the  spline  are  continuous  at  all  joints;  end  (c)  the 
lifted  to  satisfy  new  pr  additional  data.  The  experience  f 
• he  waterway::  Experiment  station  (WKU)  in  applying  spline  techniques  to  engineering 
problem  la  indicated  that  these  advantages  outweigh  the  additional  storage  and/or 
o r .-np  jtH’.lon  requirements  of  the  technique  in  many  applications,  finer  the 

spline  futictii.i  r paired  to  pass  through  all  data  j:oint.s , erratic  derivative  be- 
•.av i • i :iay  re  ill.  fr  m experimental  error  when  the  data  points  are  numerous  and 
1 ely  spa  ,-d.  Trial  and  error  methyls  for  smoothing  such  functions  exist,  but  they 
nr*-  t ime  <•  i.  uming.  Wh  ' experience  ha:  indicated  tliat  acceptable  results  can  gener- 
ally be  obtain*  . by  -imply  selecting  a more  limited,  more  widely  spaced  set  of  the 
data  points  to  which  to  fit  the  curve. 
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